Stellar Dynamics of Extreme-Mass-Ratio Inspirals 
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Inspiral of compact stellar remnants into massive black holes (MBHs) is accompanied by the 
emission of gravitational waves at frequencies that are potentially detectable by space-based inter- 
ferometers. Event rates computed from statistical (Fokker-Planck, Monte-Carlo) approaches span 
a wide range due to uncertaintities about the rate coefficients. Here we present results from direct 
integration of the post-Newtonian TV-body equations of motion descrbing dense clusters of com- 
pact stars around Schwarzschild MBHs. These simulations embody an essentially exact (at the 
post-Newtonian level) treatment of the interplay between stellar dynamical relaxation, relativistic 
precession, and gravitational-wave energy loss. The rate of capture of stars by the MBH is found 
to be greatly reduced by relativistic precession, which limits the ability of torques from the stellar 
potential to change orbital angular momenta. Penetration of this "Schwarzschild barrier" does oc- 
casionally occur, resulting in capture of stars onto orbits that gradually inspiral due to gravitational 
wave emission; we discuss two mechanisms for barrier penetration and find evidence for both in the 
simulations. We derive an approximate formula for the capture rate, which predicts that captures 
would be strongly disfavored from orbits with semi-major axes below a certain value; this prediction, 
as well as the predicted rate, are verified in the A^'-body integrations. We discuss the implications 
of our results for the detection of extreme-mass-ratio inspirals from galactic nuclei with a range of 
physical properties. 



PACS numbers: Valid PACS appear here 



INTRODUCTION 



Compact stellar remnants - white dwarfs, neutron 
stars, and stellar-mass black holes (BHs) - can be cap- 
tured by massive black holes (MBHs) at the centers of 
galaxies like the Milky Way, without suffering tidal dis- 
ruption in the process. Such extreme-mass-ratio inspirals 
(EMRIs) are a potential source of low-frequency gravi- 
tational waves for space-based gravitational wave (GW) 
interferometers HHjl . Capture orbits for EMRIs can be 
very eccentric j42|, displaying extreme versions of rela- 
tivistic precession. Typical EMRIs will have low, instan- 
taneous GW amplitudes, but signals can potentially be 
observed over > 10^ cycles as the compact objects gradu- 
ally lose energy and spiral in, allowing the signal-to-noise 
ratio to be built up over time using matched filtering or 
other techniques [Sl-Q- Detailed information about the 
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structure of spacetime is encoded in the GW signal^per- 
mitting strong-field tests of theories of gravity l&--10|. 

Predictions of the EMRI event rate span a wide range, 
from ~ 10~^ yr~^ to ~ 10~^ yr~^ per galaxy 'TP ^Uf . 
There are two basic sources of uncertainty. Only stars 
originating from tightly bound orbits, a < 10 mpc (milli- 
parsecs), can complete their inspiral without being scat- 
tered prematurely into the MBH or onto a wider orbit. 
But the number and distribution of stars and stellar rem- 
nants at these radii is essentially unconstrained, even in 
the Milky Way [13, and estimates of the event rate must 
therefore be based on extrapolation of the stellar distri- 
bution observed on much larger scales, or on theoretical 
models. In addition, the collisional dynamics of relativis- 
tic star clusters around MBHs are poorly understood. 

This paper addresses the second source of uncertainty. 
We present results from long-term (10^ — 10^ yr ), direct 
A^-body simulations of clusters of compact stars around 
a MBH. Relativistic corrections to the equations of mo- 
tion are included up to 2.5 post- Newtonian order [Tsl . 
hereafter Paper I]. These new simulations permit an es- 
sentially exact treatment of the interplay between stellar 
relaxation and GW emission, avoiding the approxima- 
tions that must be made in statistical (Fokker-Planck, 
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Monte-Carlo) treatments. 

In such a statistical treatment of EMRI inspiral, Hop- 
man & Alexander [l^ have shown that the dynamical 
evolution leading to capture on an inspiral orbit is driven 
by "resonant relaxation" 19] due to the residual torques 
from the stellar background. They argued that relativis- 
tic, in-plane (Schwarzschild) precession plays a critical 
role in suppressing the stellar torques on eccentric orbits, 
thereby allowing the stars to follow quasi-periodic EMRI 
inspiral orbits, rather than be strongly torqued into di- 
rect plunge orbits; the latter would produce non-periodic, 
broad-band gravitational wave events that would be dif- 
ficult to detect. 

The simulations presented here reveal that the inter- 
play between Newtonian torques and relativistic preces- 
sion not only limits the effectiveness of stellar relaxation 
before it can drive stars into plunge orbits, but in fact 
creates a dynamical barrier, the "Schwarzschild barrier," 
which repels stars back to less eccentric orbits, thereby 
strongly mediating the EMRI rate. We develop a Hamil- 
tonian model for the behavior of orbits near this barrier 
and use it to identify two modes by which stars can cross 
the barrier and become EMRIs. Evidence for both modes 
of barrier penetration are found in the A^-body simula- 
tions. 

In Sec.|TT]we summarize the computational techniques 
and the A'^-body initial conditions. In subsequent sec- 
tions we present results from integrations in which the 
PN terms are absent, or included only at the 2.5PN (GW 
emission) level (Sec.lIVP: and in which all PN terms up to 
and including 2.5PN are included (Sec.|V]). In Sec. IVlwe 
also present an extended discussion of orbital dynamics 
near the Schwarzschild barrier based on a Hamiltonian 
formulation. Sec. I VII discusses the implications of our 
results for the rate of EMRI production in real galaxies 
and Sec. IVIII sums up. 

We confine ourselves in this paper to non-rotating, i.e. 
Schwarzschild, MBHs. The consequences of spin will be 
discussed in a subsequent paper. 

II. MODELS AND METHODS 

The A^-body systems considered here consist of a sin- 
gle massive particle, representing a massive black hole 
(MBH), and 50 lower- mass particles representing stellar 
remnants (referred to below, interchangeably, as either 
BHs or stars). Each BH particle had a mass 5 x 10~^ 
that of the MBH particle. If the latter is assigned a mass 
of 

M. = 1 X 1O^A^0, (1) 

the BH particles have masses of rn — SOA^©. This value 
is somewhat larger than the predicted masses of the BHs 
that form in stellar collapse, i.e. 10 — 20Mq [20|. The 
choice for m/M, was motivated by the need to integrate 
the A^-body systems for a time of order the two-body 
relaxation time, which scales as for a system of fixed 



N. Alternatively, if m is set to IOMq, M, = 2x lO^A^©; 
however we note that the existence of AIBHs with Af, < 
1O^A^0 is speculative. 

Unless otherwise stated, we adopt Af, = 1.0 x lO^A/J© 
below when quoting A^-body results in physical units. In 
most cases, the dynamical theory used to interpret the 
A^-body results will allow the event rates derived here to 
be scaled approximately to systems of different m and 
M.. 

The initial orbital elements of the BH particles were 
selected randomly from the distribution 

N{a, e'^)dade^ = Nodade^ (2) 

with a and e the semi-major axis and eccentricity of the 
Keplerian orbit about the AIBH. Eq. 1^ corresponds to 
an isotropic (in velocity) distribution with configuration- 
space density 

n{r) = no ( — ) . (3) 

This is roughly the expected radial de pen dence for a re- 
laxed population around a MBH [H [HI, [13] . The ini- 
tial distribution in semi-major axis was truncated above 
a = 02 = 10 mpc and below a — ai =0.1 mpc. Setting 
A^ = 50 and m = 5OA^0, the enclosed, distributed mass 
becomes 

M4< r) w 25OA^0f , f < 10 (4) 

where f is the radius in units of mpc and the subscript 
"*" indicates the distributed mass, i.e., the stars. 

While the values of A^ and m /AI, were chosen primar- 
ily on the basis of computational convenience, the models 
adopted here are not necessarily poor representations of 
real galactic nuclei. Steady-state models of the center 
of the Milky Way galaxy [jj, fl(| typically find that the 
distributed mass within ~ 100 mpc of the MBH is dom- 
inated by stellar BHs (as opposed other types of stellar 
remnant, or stars) with M^,(r < lOmpc) w lO'^A/f©. Ex- 
pressed in terms of the gravitational radius defined in 
Eq. pI7|) . and assuming p oc , a distributed mass of 
IO'^Mq within 10 mpc implies 

M.(<.)« 200 (^). (5) 

By comparison, the scaling adopted above implies that 
for our models, 

M.i< r) « 120 (^) . (6) 

The A^-body integrator is described elsewhere [1^, [13] ; 
it includes post-Newtonian accelerations of orders up to 
and including 2.5 (i.e. c~^) in the interactions between 
the MBH- and star particles. The algorithm was modi- 
fied for the current study to allow merger of star particles 
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with the MBH. The condition for a merger was an instan- 
taneous separation [25] 



r < Tcapt = 8rg 



(7) 



or ~ 4 X 10"'' mpc if M, = lO^A^©. The angular mo- 
mentum and eccentricity of an orbit that just grazes the 
capture sphere are (in the Keplerian approximation) 



capt 1 



1 



^capt 



(8) 



For ai < a < a2 and using the adopted value of rcapt, 
the eccentricity of a capture orbit satisfies 



4 X 10"^ < 1 - < 4 X 10 



-3 



(9) 



The mass of a merged particle was added to that of the 
MBH in such a way that linear momentum was con- 
served. 

An EMRI event was defined as any merger occurring 
from an orbit with semi-major axis, at the moment of 
capture, less than 0.01 mpc. Mergers occurring from or- 
bits with larger a were recorded as "plunges" . 

While capture in its final stages would be driven by 
energy loss due to emission of gravitational waves, as 
represented here by the 2.5PN terms, the capture rate is 
determined primarily by dynamical interactions that take 
place far beyond rg. In order to better understand the 
dynamical mechanisms leading to capture, three series of 
iV-body integrations were carried out, incorporating dif- 
ferent subsets of the full PN equations of motion defined 
in Paper I. 

Series I: No PN terms were included. Stars were nev- 
ertheless allowed to merge with the MBH if they passed 
within rcapt ("plunges"). 

Series II: The 2.5PN terms were included. As a re- 
sult, some stars ("EMRIs") were captured onto orbits for 
which the timescale for gravitational wave energy loss is 
less than the timescale for scattering by other stars. 

Series III: AU PN terms (IPN, 2PN, 2.5PN) were in- 
cluded. 

In each series, at least eight different Monte Carlo re- 
alizations of the same initial conditions were integrated 
forward in time. Models from Series I and II were inte- 
grated for a time of 10^ yr, based on the scalings adopted 
above. For models from Series HI, inclusion of the addi- 
tional PN terms caused the integrator to run more slowly, 
and most integrations were terminated after 2 — 3 Myr. 
Calculations were carried out on gravitySimulator, the 
32-node cluster at RIT. 



III. BASIC SCALES OF LENGTH AND TIME 

Here we define length and time scales associated with 
an idealized model consisting of a central MBH and 



a smooth, spherical distribution of surrounding stars. 
Other time scales, associated with coUisional (relaxation) 
processes, are defined below. 

The length scale associated with the event horizon of 
the MBH is 



GM, 



4.80 X 10^^ 



mpc. 



(10) 



where the numerical value assumes M, ~ 1.0 x lO^A^©. 

Ignoring the contribution of the stellar BHs to the 
gravitational potential, the (Newtonian) orbital period 
of a test mass around the MBH is 



Pr 



2^3/2 



vgm: 



2.96 a^^^yi 



(11) 



where a is the test mass's semi-major axis in units of mpc; 
the second relation again assumes M, — 1.0 x IO^/Aq. 

Approximating the stellar BHs as a smooth mass dis- 
tribution, p{r) = mn{r) with n{r) given by Eq. ([3]), their 
contribution to the gravitational potential is 



$*(r) = In — 

ro \r() 



constant 



(12) 



where Mi,(< r) — Mo(r/ro); setting tq — 1 mpc gives 
Mq — 2507^0 in our models. 

Deviation of the potential from that of a point mass 
induces a precession in the (fixed) plane of an orbiting 
star, the "mass precession." Orbital perturbation theory 
[e.g. gives for the precession rate, in the limit Mq <C 
M., 



duj 
H 



M^{r < a) VT^e^ 



M, 



1 + 



(13) 



where 



— = (GAf.)i/2a-3/2 



(0.47yr)-^5 



(14) 



is the radial frequency and e is the orbital eccentricity. 
Precession is retrograde, i.e. in the opposite sense to the 
orbital motion. For the adopted mass model [43| . 



5(e) 



1 + Vl~^ 



2VT 



1.18 X lOV) 5(e)a^/^ (15a) 

(15b) 



In the limit e — > 1, Eq. predicts vm 0, i.e. radial 
orbits do not precess. 

The post-Newtonian accelerations also contribute 
to the in-plane precession. To lowest order, the 
Schwarzschild contribution is 



3 {GM, 



,3/2 



3r„ 



c2 (1 - e2) 0^5/2 



a(l-e2) 



(16) 
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in the opposite (prograde) sense, and 



^GR 



(1.02 X lOV) (1 - e^) 



2^ 5^/2. 



(17) 



While we defer a detailed treatment of spin effects to a 
subsequent paper, we note here that the Kerr contribu- 
tion to the in-plane precession is smaller than ([TBI) by 
a factor ~ x{'''g/p)^^'^^ where x is the dimensionless spin 
parameter of the MBH and p = {\ — e'^)a is the semi-latus 
rectum. Excepting very shortly before a merger, this fac- 
tor would be much smaller than unity in our simulations. 



IV. SERIES I AND II 

As discussed in more detail below, including the rela- 
tivistic terms in the A'^-body equations of motion resulted 
in much lower EMRI rates than expected based on New- 
tonian dynamics of a compact cluster around a MBH. 
The essential element that differs between the relativis- 
tic and non-relativistic dynamics turns out to be the IPN 
precession of the periapse, Eq. (ITBl) . In order to quantify 
the magnitude of the differences, two sets of experiments 
were carried out in which some or all of the relativis- 
tic terms were omitted. Integrations from Series I were 
based on the Newtonian equations of motion. Series II 
included also the 2.5PN terms, allowing capture of stars 
onto inspiral orbits via GW energy loss. In integrations 
from both Series I and Series II, the IPN (Schwarzschild) 
precession is absent. 



In this (Newtonian) regime, the mechanism expected 
to dominate the scattering of stars onto high-eccentricity 
orbits around a point mass is resonant relaxation (RR) 
[igj . Because the smooth gravitational potential has 
symmetries that restrict the orbital evolution (i.e. to 
fixed ellipses in the case of a Newtonian point mass), 
perturbations on a test star are not random but corre- 
lated. This leads to coherent changes AL — Tt on times 
i ^ icoh by the residual torque |T| ft! ^fNGm/r exerted 
by the N randomly oriented, orbit-averaged mass distri- 
butions of the surrounding stars. The coherence time is 
set by the mechanism that most rapidly causes orbits to 
precess, randomizing T. In these non-relativistic simula- 
tions, that mechanism is mass precession, Eq. ()15a|) . The 
accumulated change over icoh, |ALcoh| ^ |Ttcoh|, then 
becomes the mean free path in L space for the long-term 
(t ^ tcoh) random walk in L. The effective relaxation 
time associated with RR satisfies 



lAL, 



coh| 



t 

tcoh 



I.e. 



tRR — 



V AL, 



coll 



t 

tRK 



t-coh 



(18) 



(19) 



where Lc = V GAf, a is the angular momentum of a cir- 
cular orbit with radius r « a. These relations should 
be understood as correct only in an order-of-magnitude 
sense. 

In the coherent regime the change in orbital angular 
momentum is 



A. Series I 

Fig. [1] shows, for one integration in Series I, the evo- 
lution of semi-major axis a and eccentricity e for each of 
the 50 stars, until a time of 2 Myr. Star-star gravitational 
scattering induces substantial changes in the orbital an- 
gular momenta over these time scales, while the energy 
(i.e. semi-major axis) remains nearly constant. 

Whenever the periapse distance Vp = a(l — e) drops 
below Tcapt the star is captured. Almost all such events 
are "plunges" since no GW energy loss occurs, and since 
essentially no stars have initial semi-major axes less than 
0.01 mpc, the condition defined above for a capture to be 
classified as an EMRI. In the simulation shown in Fig. [TJ 
17 out of the 50 stars are captured by t ~ 2 Myr and 30 
stars are captured by t = 10 Myr. 

Fig. [5] shows the time-averaged capture rate as a func- 
tion of time, defined as the number of mergers occurring 
in time t divided by t. Events from each simulation in 
the series were summed and the result was divided by the 
number of simulations; in other words, the plotted rates 
refer to a cluster with N{t = 0) = 50. The capture rate 
drops with time, since both the number of stars available 
to merge, and the number of stars acting as scatterers, 
decrease with time. 



dL 



dt 



N m 



I — Gm Lr 



(20) 



where N is roughly the number of stars within a sphere 
of radius a, and /3s is a dimensionless factor of order unity 

mill. 

The precession time due to the distributed mass, 
Eq. (|15a|) . can be written 



tu ~ (1.2 X 10V)5(e)ai/' 



/ M. \ / Af, 



V 106X0/ \250Mq 



(21) 

where Mi,{< r) = Mof and Mq = 25OA^0 for the models 
considered here. 

Identifying with tcoh, and writing N{< r) — A/^(< 
r)/m, Eqs. give 



M 
ni 



5.9 X 10-* 



yr / M, 



1/2 



/ TO 

Uox© 



&3 



(22a) 
(l2b) 



We note that has dropped out. This is only valid for 
values of M^, large enough that tcoh — tjM is shorter than 
all other time scales of interest. 
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FIG. 1: A simulation from Series I (Newtonian), a and e are the semi-major axis and eccentricity respectively of the two-body 
system consisting of one star and the MBH; rpcri = (1 — e)a and Tapo = (1 -|- e)a. Different colors correspond to different 
particles (since the total number of colors available was 12, each color is used for more than one particle). Dashed lines indicate 
the capture radius, rcapt = 8rg. 



In the expression for iRR, the form of the density pro- 
file is still reflected in the dependence of g on e, which, 
for the assumed initial mass distribution, is strongly de- 
pendent on e as e — >■ 1 (Eg. I15b|) . However what matters 
for the coherence breaking is the relative precession of 
the test particle's orbit with respect to the other orbits. 
Hence it is reasonable to average g{e) in Eq. ([^ over 
the eccentricity distribution for all the stars, Eq. ([2]): 

J\{e)de'^^. (23) 



In what follows we ignore changes in g due to evolution 
of M4r). 



Because the integration time is long compared with 
tRR, we expect a quasi-steady-state to be set up in the 
angular momentum distribution at each a, such that the 
rate of loss of stars into the capture sphere is roughly 
equal to the rate at which new stars are being scattered 
onto low-angular-momentum orbits. In this regime, the 
(differential) rate at which stars are scattered into the 
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time (yr) 



FIG. 2: Time-averaged capture rates, defined as the total number of events until time t divided by t, computed from the 
complete set of runs in each series. Series I: dashed line is the prediction of Eqs. (|26 p -(|27 p for Ci = 0.5. The solid (black) 
curves in the other panels are the total event rates, EMRI plus plunges. 
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FIG. 3: Open circles show the number of stars captured up 
until time t in the combined runs from Series I; values are 
normalized to an initial total number of 50. Plotted events 
are all "plunges." Solid lines show the predictions of Eqs. (|26|l - 
([27ll for Ci = (0.3,0.5,0.7). 



MBH at each a is given approximately by: 

N{a,t)da 



r(a, t)da ; 



(24) 



[e.g. 12911. The logarithmic term can be interpreted as 
the approximate number of relaxation times required for 
an orbit to diffuse in angular momentum from e to 
e « e,n [aOl. 

Using Eq. (P^. the differential capture rate can be 
written 



r(a, t)da — Ci 



N{a,t)da 



M, In (1 



=2 ^-l/2 



Pia) 



(25) 



where Pfg and all other uncertainties have been absorbed 
into the fitting parameter Ci , assumed independent of a 
and t. 

Equating T{a,t) with —dN{a,t)/dt and using Eqs. ^ 
and (|lip . the evolution equation for N{a,t) becomes 



ON 
'd7 



iV(a,T) 



a3/2ln(a/rcapt)' 



T = t/ti, (26a) 



h = (5.9 X lOVOCi" 



1 M,/m / M. 

2 X 104 [WM^ 



-1/2 



(26b) 



with solution 



N{d,T) = 7V(5,0)e~^/^S 
Ti = a^/^ In (a/rcapt) . 



(27a) 
(27b) 



Fig. 13] plots the predicted, cumulative number of 
events versus time, compared with the results from 
the Series I integrations. The agreement is good 
for Ci ~ 0.5. The mean capture rate is given by 
i-i /a7 [^("' 0) ~ ^("' 0] this is plotted, with Ci = 
0.5, as the dashed line in Fig. [2^. 



B. Series II 

Series 11 integrations included the 2.5 PN terms in the 
equations of motion, allowing some stars to be captured 
onto orbits that inspiral gradually into the MBH via GW 
energy loss. 

The energy loss timescale associated with the 2.5PN 
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-0.5 0.5 

log [a(t=0)/mpc] 



FIG. 4: Distribution of initial semi-major axes for the cap- 
ture events from Series II. Red (unfilled) histogram shows the 
plunges; blue (cross-hatched) histogram shows the EMRIs. 
Elapsed time is 10*^ yr. 



with the RR timescale defined above. Equating 
with icw then gives the condition for capture onto an 
inspiral orbit: 



a(l 



1 /3407r\ 

2 ['^J 



2/5 



(32) 



i.e. capture requires a periapse distance of ~ Sr^. This 
is shghtly smaller than the separation at which mergers 
were assumed to take place in the simulations (Eq. [7|) . 

Given the approximate nature of Eq. p2p . one expects 
capture onto inspiral orbits for some fraction, of order 
unity, of stars that would otherwise plunge into the MBH. 
Fig. U shows the cumulative histogram of capture events 
for the integrations of Series II. Roughly 1/4 (54 out of 
206) events were EMRIs. Time-averaged capture rates 
are shown in Fig. ^p. These results are consistent with 
expectations. 

As shown in the next section, these results are sub- 
stantially changed by the inclusion of the IPN and 2PN 
terms into the equations of motion. 



SERIES III 



terms is 13 ll 



1 da 
a dt 

5 ( 



64 G^Ml 



(1 



,7/2 



73 



37 



1 + + 



1.2 X 10 



14 



507W. 



24 
M, 



96 

-2 



(1 



where the latter expression assumes e ~ 1. In this limit, 
GW inspiral occurs along lines of fixed slope in the a, (1 — 
e) plane until shortly before the merger: 



A(l 



Aa 



1 - e 



(29) 



such that Tp = (1 — e)a is approximately constant [3l|- In 
order to avoid plunging, a star must reach a high enough 
eccentricity that the GW timescale is shorter than the 
time for gravitational encounters to scatter the star onto 
a different orbit. 

From Eqs. pVB|) -(l^. the time required for GWs to 
change e by of order 1 — e is ^gw- In the case of 
gravitational encounters, changes in angular momentum 
are equivalent to changes in eccentricity since a is nearly 
conserved. The time for encounters to change L by of 
order itself is 



Fig. [S] shows the evolution on the (a^, 1 — e^) plane of 
all 50 stars in an integration from Series III. Integrations 
(28)n this series included all PN terms (IPN, 2PN, 2.5PN). 
The quantities the IPN generalizations of the 

Keplerian semi-major axis and eccentricity respectively; 
to this PN order, the periapse and apoapse distances are 
given respectively by (1 — 6^)0^ and (1 -I- er)ar [32j . 
g-j7/2 Over the course of the 2 Myr interval plotted in Fig. [51 
only one capture occurs: an EMRI. The plunge events 
that dominated the integrations from Series I and II are 
absent. The mean capture rates computed from all in- 
tegrations in Series III are shown in Fig. [51 The mean 
capture rate is < 1 Myr~^, with 74% of the events EM- 
RIs. By comparison, in Series I and II, the mean capture 
rate was > 10 Myr~^ at early times, and almost all events 
were plunges. 

The proximate reason for the much lower event rate in 
the integrations from Series III is suggested by Fig. [SI 
there is an apparent barrier in orbital eccentricity, or 
angular momentum, which very few stars cross on Myr 
timescales. Furthermore, the single star that is captured 
in that figure appears to require a time much longer than 
~ (1 — er)iRR to cross the gap. The barrier is illustrated 
more clearly in Fig. [6l based on another integration from 
Series III. This figure shows that there is some "barrier 
penetration" for orbits with small semi-major axis; we 
discuss the reasons below. 



(AL)^ 



tRR 



(30) 



The Schwarzschild barrier 



L 

-J— I ^RR 



2(1 - e)iRR 



(31) 



Adding the IPN terms to the equations of motion re- 
sults in precession of the argument uj of orbital periapse. 
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5x10^ 10^ 1.5x10^ 2x10^ 

time (yr) 

FIG. 5: A simulation from Series III (all PN terms), ar and are the IPN generalizations of the semi-major axis and 
eccentricity; rperi = (1 — er)a,- and rapo = (1 + er)ar. Different colors correspond to different particles (the number of different 
colors is 12 so each color is used for more than one particle). Dashed (black) lines show the assumed capture radius, rcapt ~ 8rg. 
In the top frame, the dotted (red) line is the Schwarzschild barrier, Eq. (|36|l , and the dash-dotted (blue) line is the approximate 
condition for GW capture, Eq. 



with an orbit-averaged frequency given by Eq. (jl6|) . For 
low eccentricity orbits, the rate of this Schwarzschild pre- 
cession Q is comparable to that produced by the dis- 
tributed mass, Eq. (IT3t . at the radii of interest here. 
But whereas the latter rate tends to zero as e — > 1, the 
Schwarzschild precession rate diverges, as ~ (1 — e)~^. 
The effective time over which background torques can act 
is determined by the fastest mechanism that changes the 
relative orientation of a star with respect to the gravita- 
tional field produced by all the other stars. For a highly 



eccentric orbit, this mechanism is Schwarzschild preces- 
sion and its associated timescale tends to zero as e — ^ 1. 

We suggest that the angular momentum barrier be 
identified, in a qualitative way, with the value of L at 
which the torques become ineffective due to the orbit's 
rapid Schwarzschild precession. 

The residual torque produced by an otherwise- 
spherical distribution of stars, at r « a, is of order 



Gm 



1 



GM,(a) 



(33) 
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FIG. 6: Illustrating the angular momentum barrier when all PN terms are included, in two short time segments extracted from 
a Series III integration, (a) 4 x lO^yr < t < 8 x 10^ yr; (b) 1 x lO'' yr < t < 1.4 x lO'' yr. Dashed (black) line is the capture 
radius; dotted (red) line is the predicted angular momentum barrier, Eq. (|36|) . Stars that lay initially to the left of the barrier 
were excluded. This integration produced no EMRIs. The time interval plotted. At = 4 x 10^ yr, is somewhat longer than the 
RR timescale of Eq. (|22p and much longer than (1 — er)tRR. Note the "barrier penetration" at small values of Ur- 



where Mi,{a) = mN(a) is the distributed mass within 

1/2 

radius r = a. Writing L — [GM,a{l — e^)] for the 
angular momentum of a test star, the time scale over 
which this fixed torque changes L is 



1 dL 



M, 

M{a) 



a3(l 



GM. 



1/2 



(34) 



The condition that this time be shorter than the rela- 
tivistic precession time, it /vq^, is 



M, 



SB 



(35) 



where we have written (. = L/Lc = (1 — e^)^/^. Evalu- 
ating the quantities in Eq. ([55)1 for the A^-body models, 
the critical semi-major axis becomes 



CsB (1 - e') 



-1/3 



(36) 



where Csb is a constant of order unity. Eq. ([36]), with 
CsB = 0-7, is plotted as the dotted (red) lines in Figs. [SI 
[l[7]and[ll 

Assuming that the condition (l36l) holds for all values 
of a, the normalizing constant Csb can be interpreted as 
the value of a when e = 0, i.e. as the minimum value of 
a for which the barrier exists. One expects that orbits 
with semi-major axes larger than this minimum value 
(and smaller than a maximum value, defined below), and 
that approach the barrier from the right on Figs. [5] or [51 
will have a hard time crossing it, since torques become 
inefficient near the barrier. 

We develop these ideas more quantitatively in the next 
subsection. Before doing so we present a more quantita- 
tive model for the behavior of low-angular-momentum 



orbits under the combined influence of relativistic pre- 
cession and Newtonian torques. 

Fig. [7^ shows the evolutionary track of a star from a 
Series III integration. The star first strikes the barrier at 
t « 1.8 x 10^ yr; the eccentricity then oscillates several 
times at roughly fixed amplitude before decreasing again, 
carrying the star away from the barrier at t > 2.2 x 10^ yr. 
During each bounce, the argument of periapse lo advances 
by ~ 27r. 

Many other examples of "bounce" near the angular 
momentum barrier were extracted from the TV-body in- 
tegrations. While differing in detail, all such orbits ex- 
hibited a variation in eccentricity near the barrier with 
a period roughly equal to the period of Schwarzschild 
precession. 

This feature suggests that the torques responsible for 
angular momentum changes near the barrier are due to a 
distortion of the stellar potential, expressed about the lo- 
cation of the MBH particle, that is lop-sided or dipole in 
character. (A quadrupole distortion would cause changes 
in L at twice the frequency of circulation, etc.). That 
the dominant component of the torquing potential should 
have such a form is not unreasonable, since if one repre- 
sents the gravitational potential from TV orbit-averaged 
stars in a multipole expansion, the largest terms are ex- 
pected be the monopole (due to the spherical cluster) 
followed by the dipole (due to y/N departures from spher- 
ical symmetry) etc. 

We also verified that the behavior of orbits like that 
plotted in Fig. [7] was unchanged if the mass of the test 
particle was drastically reduced. These tests confirmed 
that the variations in orbital angular momentum in the 
test particle's orbit were not a spurious result of motion 
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of the black hole particle induced by the test star's pre- 
cession. 

We used the following simple potential to model the 
motion of a test star subject to a lopsided force from all 
the other stars: 

^(r) = -^^ + ^,Vir)- S{a)acose. (37) 
r 



The second term on the right hand side of (|37p is the 
potential of the spherical star cluster. For the models 
considered here, 



V{r) = In 



ro 



ro 



Mo ^M^r < ro). 



(38) 

The third term represents the lopsided distortion of the 
stellar potential; the amplitude S*, which has dimensions 
of acceleration, is assumed to depend on the test-orbit's 
semi-major axis as 



S{a) 



Gmy^N(a) GM^{a) 



The corresponding density is 



PD{r,e) 



S{a)a cos 9 
27rG r2 ' 



(39) 



(40) 



The integrated mass corresponding to the lopsided com- 
ponent is zero. 

Due to the dominance of the first term in Eq. (|37| , the 
radial period of an orbit is shorter than all other orbital 
time scales. In Appendix El we express the Hamiltonian 
corresponding to the potential ([57)1 in terms of Delaunay 
action-angle variables and average over the the radial mo- 
tion, including the orbit-averaged term that generates the 
Schwarzschild precession. The result is a set of four equa- 
tions that describe the rates of change of the (osculating) 
Keplerian elements (L,L2,w,r2): 



du) 



r' - AMe/{i + e) 

A]j sin a; 



£ . 

sm i 

e 



sin i P 



— — = —Anesinicosuj, 



dr 

dn 

(M 
~dT 



— = -Ane 



^ = 



£z sinw 
P sini ' 



(41a) 
(41b) 
(41c) 
(41d) 



where we have defined the dimensionless elements i = 
L/I — Vl — e^, £z = Lz/I and cosi = £z/£ in terms of 
the radial action / = (GAf.a)^/^, and the dimensionless 
time is T = v^t^ where 



vo = Vr- 



ZGM, 



(42) 



In defining the Delaunay variables, the plane of refer- 
ence has been taken to be the (x, y) plane and the ref- 
erence direction is the x-axis; thus an orbit in the (x, z) 
plane has sinz — 1, and for such an orbit, w = 7r/2 
corresponds to an orientation parallel to the 2;-axis, the 
assumed direction of the lopsided distortion. In the case 
of an orbit in the {x^y) plane, sini = 0, £z = £■, and the 
orientation of the orbit is determined by w-l-fi; according 
to Eqs. (|41ap and (|41cl) . the terms in {d/dT){uj + ft) that 
are proportional to A]j sum to zero in this case. 

The dimensionless parameters Am and Ajj specify the 
strength of the spherical and lopsided components of the 
distributed mass: 



Am 
Ad 



1 Mi,{a) a 
3 
1 



M. 

S 



1 M^{a) a 



(43a) 



(43b) 



We note that Am/Ajj sa N^^'^ which is of order unity in 
the models considered here. Thus whenever it is relevant 
to neglect Am, Ad is also negligible. 

After averaging, the first and third terms in Eq. ((57)) 
result in the same equations of motion as in the classical 
Stark problem [e.g.|33|. The corresponding solutions [e.g. 
[3^ consist of circulation in £7, with period 



4tt 

Stark — 777^ 

3d V a 



GM, 



(44) 



oscillations in i, e and a; have the same period, while 
the 2-component of the angular momentum is fixed. The 
eccentricity reaches a maximum value that depends on 
Lz; for Lz = 0, i.e. for i = it/2, emax = 1, while for 
Lz the maximum eccentricity is less than one. 

In the case considered here, precession of a sufficiently 
eccentric orbit is dominated by the Schwarzschild term, 
the first term on the right hand side of Eq. (|41ap . Such an 
orbit circulates in a nearly fixed plane and the eccentric- 
ity varies with a period equal to the period of circulation. 
If £ is sufficiently small, we can write duj/dr « {£)^'^ with 
{£) the time-averaged angular momentum. Eqs. (|4T|) then 
have approximate solution 



1 - 



(£)Ad sin i cos (vt) , 
3 {GM.f^^ 



(45a) 
(45b) 



In this limit, the amplitude of the angular momentum 
oscillations, 



2{£fAu sini, 



(46) 



decreases quadratically with {£) . 

We now return to the full equations of motion (|4T|) in 
order to test whether the detailed behavior of stellar tra- 
jectories near bounce in the A'^-body integrations is con- 
sistent with our simple model. In the A^-body models. 
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FIG. 7: (a) A "bounce" orbit extracted from a Series III A''- 
body integration. Plotted are the semi-major axis, argument 
of periapse and eccentricity versus time, at low (upper) and 
high (middle) time resolutions. In the plots of eccentricity 
vs. time, the lower (red) curves show the predicted location 
of the angular momentum barrier, Eq. (|36|) . with Csb = 0.7. 
Changes in the predicted barrier location reflect changes in 
the semi-major axis, (b) A solution to the equations of mo- 
tion (|4ip that reproduces the important features of the A'- 
body orbit in (a). The duration of the "bounce" phase is 
roughly the coherence time for the background (stellar) po- 
tential. Additional details are given in the text. 

the dimcnsionless parameters that appear in the equa- 
tions of motion are 

1.8&2, (47a) 
1.2&3/^ (47b) 
(3.26 X 10V)"^a"^^^ 



Am 
Ad 

^0 



(47c) 



and the Schwarzschild precession period is 



Pgk = 



2.1 X lOV (1 - e^)a^^^- (48) 



The A^-body orbit in Fig. [7^ exhibits ^ 6 full circulations 
in w in a time ~ 3.5 x 10'* yr, corresponding to a preces- 
sional period of ^ 6 x 10^ yr. The semi-major axis for 
this star during the bounce is 3.5 < a < 5 and the eccen- 
tricity is —2.8 < log2o(l — e) < —2.1. Inserting a = 4 and 
logio(l - e) = -2.5 into Eq. dM]) gives Pgr « 4 x 10^ yr 
which is quite consistent with the observed precessional 
period. 

Fig. [71d shows a solution to Eqs. (HTI) that reproduces 
the other important features of the iV-body orbit near 



0.01 



% 0.5 



0.2 





1 1 1 1 


1 1 ^J. 1 1 1 _ 


- 


^crit 

o 

o 

o 


.•' L 

r 

-♦- -. 

r 

r. 
". 

. 


o 

o 

o 


c 


^min • 


o 

o 

o o 
o 

° o 
o o 
o ^ ° 
O 

o o 
> o 


o o o 

i-H \ — 


- 

1 
I 

■ 

1 1 — 1 1 1 1 


-e 1 1 1 1 1— 1- 

'"""^OOOo n 

^ O 

o 


1 

o 


\ 




, , 1 





0.01 



0.1 



FIG. 8: Properties of two-dimensional (sini = 7r/2) solutions 
to the equations of motion (|4ip , with Am ~ 30, Ao = 10. {£} 
is the dimensionless angular momentum, £ — \/l — e'^, aver- 
aged over one precessional or librational period P. Top panel: 
minimum and maximum angular momenta reached by the or- 
bit over one period. Filled circles ("L") are librating orbits 
while open circles ("C") are circulating orbits. The dashed 
lines at ^ = iait mark the angular momentum at which the 
precession rate, duj/dt, is zero. The dotted line marked l-mi-n is 
the minimum angular momentum reached by librating orbits; 
it is argued in the text that this is essentially the angular mo- 
mentum corresponding to the Schwarzschild barrier, Eq. p5|) . 
Bottom panel: Librational/precessional periods as a fraction 
of the Schwarzschild period, computed from Eq. [16] after re- 
placing (1 - e^) by (£)^ 



bounce. We set A-q = 10 and A-q — 30 i.e., a « 4; the 
initial values of the orbital elements were log(l — e) = 
—2.1, Lu = — 7r/2, = 0, and i — 0.357r. Variations in 
and i (not shown here) were similar in amplitude for the 
A^-body and numerically computed orbits. 

We carried out similar comparisons for other orbits 
near the barrier. While differing in details, all the cases 
examined could be adequately represented via solutions 
to our simple Hamiltonian (jA6l) . 

The assumptions made in deriving the Hamilto- 
nian (jA6[) are not specific to orbits near the barrier. Fig. [8] 
summarizes the properties of orbits, of arbitrary angular 
momentum but restricted to the x — z plane (sin i = 7r/2), 
in the potential of Fig. [7}d. Orbits can either circulate 
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FIG. 9: Two orbits from FigM (a) {^) = 0.225 (circulating); 
(b) {£) — 0.402 (librating). cj = n/2 corresponds to orien- 
tation along the z axis. The dotted lines in the right-hand 
panels indicate imin- 



(small angular momentum) or librate (large angular mo- 
mentum). There is a critical value of the angular mo- 
mentum, at the time the orbit is oriented parallel to the 
z-axis (i.e. the direction of the lopsided distortion), such 
that the precession rate ij — and the orbit remains fixed 
in orientation; from Eq. ()41a[) this occurs when £ = .^crit 
where 



= 1-^ 



M 



£^ ■ 

^crit 
1 + ^crit 



At 



''cnt 



^crit 



(49) 



or ^crit ~ 0.310 in Fig. |S1 Away from this value, librat- 
ing orbits experience both their minimum and maximum 
angular momenta when precessing past the z-axis; first 
in one sense (when the mass precession term dominates) 
and then in the other (when the Schwarzschild term dom- 
inates). Libration changes to circulation when the orbit 
precesses by an angle ±7r from its starting point along 
the z-axis. The minimum angular momentum reached 
by the orbit at the libration/circulation boundary is la- 
belled £min on Fig.[S] in this potential, ^min ~ 0.045. 

Fig.[S]shows two orbits from this plane, one circulating 
and one librating. We note here one property common to 
both types of orbit: stars tend to spend more time with 
high angular momentum than with low angular momen- 
tum. The reasons for this behavior are apparent in the 
equations of motion (HT|) . (1) For large £, M/dr is small, 
i.e. the orbit-averaged effects of the torque are small. (2) 
Orbits tend to linger at values of lo corresponding to large 
(. since the Schwarzschild precession rate is proportional 
to £-2. (The latter trend reverses for orbits so circular 



that mass precession dominates the Schwarzschild pre- 
cession.) 

Combining Eqs. ^ and (|43bD . 



€sB ~ (2Ad)-^ 



(50) 



or £sB ~ 0.05 for = 10. Not coincidentally, this 
is roughly equal to £min ~ 0.045. Fig. [S] shows that 
^min specifies not only the minimum angular momen- 
tum achievable by librating orbits, but is also roughly 
the minimum i reached by orbits whose angular momen- 
tum changes by of order itself over one period; these were 
the two assumptions made in deriving (I35p . 

At the same time, it is clear from Fig. [5] that orbits 
with £^ <^ (mill do exist. Apparently, such orbits are 
rarely reached in the A^-body integrations. We discuss 
the reasons in the next sub-section. 



B. Barrier penetration. I. 

Here we address the question of how orbits evolve after 
striking the angular momentum barrier. The mechanism 
( "tunneling" ) explored in this section, which is based on 
resonant relaxation, will turn out to be less important 
as a source of barrier penetration than the mechanism 
presented in the next section, based on non-resonant re- 
laxation. We nevertheless explore it in some detail since 
doing so will lead to insights about why the barrier is so 
"hard" on time scales comparable to the RR time. 

The Hamiltonian model just presented assumed a fixed 
gravitational potential. The (constant) term responsi- 
ble for the torques was assumed to arise from the timc- 
averaged potential of the N stars. In reality, the back- 
ground potential must change as the orbits of all the stars 
evolve (e.g. precess), leading to quasi-random changes in 
the direction and amplitude of the torque that acts on a 
single star. The changes in the background torques are 
responsible both for moving a star toward the barrier and 
moving it away. As we show here, such changes can also 
result in barrier penetration. 

In Sec. llVi the mass precession time. 



250Mc 



-1 

(5ila) 



was taken as the "coherence time" over which the back- 
ground torques can be assumed constant. We empha- 
size that the relevant time here is the precession time 
for a typical orbit, hence we have set g{e) ='g — 3/2 in 
Eq. (|21[) : the fact that some (high-e) orbits precess much 
faster due to relativity is not important. 

A second relevant time scale is the so-called vector res- 
onant relaxation (VRR) time, the time for orbital planes 
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to change due to their mutual torques [e.g. [lj| : 

TT M. 

tRR.v 



1.3 X lOV 




where we have written N — N{< a) /a; N ^ 5 in the 
iV-body models. Since ~ iRR,v/'\/^, one normally 
assumes ^ %r,v However the small values of N 
considered here, together with the approximate nature 
of Eq. (15^ . means that the two time scales are essen- 
tially the same in these models at all radii of interest. 
Furthermore, VRR leads to full randomization of orbital 
orientations in the sense that it changes orbital planes 
as well as orientations within the plane; mass precession 
leaves the orbital planes unchanged. For these reasons, 
we adopt ^rr^v as the coherence time in the remainder of 
this section. 

The time-independent model presented above implic- 
itly assumed that fcoh was long compared with orbital 
precessional periods. In fact, the ratio of the coherence 
time to the Schwarzschild precession time is 



^coh 
*GR 



r„M. 1 1 



a m 



ivi 



Using Eq. (|35|) to relate a to e along the barrier. 



^coh \ 
iCR / SB 



a m I — 
3— — V7V 



7.2a3/2 



(53a) 
(53b) 

(54a) 
(54b) 



where the final expression refers to the A^-body models 
and uses Eq. (j47b|) . Eqs. ([54| suggest that for orbits 
near the barrier with a > 1, the background potential 
should remain constant for several Schwarzschild preces- 
sional periods. This is consistent with the observed be- 
havior of orbits like the one in Fig. [7^ (a « 4), which 
precesses a few times before (presumably) changes in the 
background potential cause the orbit to evolve away from 
the barrier. In the case of orbits with a < 1 the two time 
scales can be assumed to be comparable in our models. 

One way to penetrate the barrier is suggested by Fig.|8l 
Random changes in the background potential, e.g. in the 
direction of the lopsided term, could have the effect of 
moving orbits progressively to the left and downward on 
that plot. This is because the minimum angular mo- 
mentum reached by an orbit over a precessional period, 
depends both on the instantaneous value of ^, and 
on the relative orientation of the orbit and the torquing 
term. A sequence of correlated changes in the direction 
of the torque could result in gradual transition down the 
narrow "neck" at the lower left of the diagram, toward 
arbitrarily small values of {l). 



This mechanism can be simply modeled if we assume 
that (1) changes in the background potential are instan- 
taneous, separated by time ^ tcoh, and (2) tcoh ^ ^gRi so 
that the orbital phase is essentially random at the time 
that the potential changes. The first assumption is not 
likely to be satisfied in all cases and we relax it below; 
however we will argue that it corresponds to the highest 
probability for barrier penetration. 

Consider first the two-dimensional case, i.e. sini — 
-K 12. Assume as well that the orbit is sufficiently far down 
the "neck" that the angular momentum follows Eq. (|45p , 

("{uj) « (^)O [1 - (£)"^D cos(w - 



(55) 



where is the initial orientation of the lopsided dis- 
tortion. Now let the direction of the distortion instanta- 
neously change, to w^. The new orbit follows 

(w) « ) 1 [1 - (^) lylo cos(a; - uj\,)\ . (56) 

If the change occurs when ui = uj^ then 

(€)i[l-(^)i^DCOs(a;'-WD)] = 
(^)0[l~(^)°^DCOs(a;°-c^f^)] • (57) 

The change in (^), A(£) = (lY - {e)° is 

A(^) « -(^)^Ad [cos^(l - cos A) - sin (5 sin A] (58) 



^ and A = ui^ — ajf^: in this simple 
model, both angles are random variables. Decreases in 
{£) are clearly allowed, although the amplitude of the 
step size becomes increasingly small, A{£) ^ as {£) 
decreases. 

There is an additional reason why evolution toward 
small {£) is disfavored. Eqs. ((55)1 assume a constant rate 
of circulation in w and ignore the eccentricity dependence 
of the averaged torque, Eq. (I41cp . But as noted above, 
away from the limit of small {£) , orbits violate both as- 
sumptions, and tend to linger at high values of £. As a 
result, changes in the potential are most likely to occur 
when £ is large. 

For these two reasons, we do not expect the mechanism 
discussed in this section to be effective at moving stars 
very far to the left of the Schwarzschild barrier; indeed 
we will argue in a subsequent paper [s^ that evolution 
to lower £ via this mechanism is exponentially suppressed 
(and this is the basis for assigning the name "tunneling" .) 
However, the ineffectiveness of RR at breaching the bar- 
rier is important in explaining why the barrier is observed 
to be so "hard," at least on time scales comparable to 

We tested this model of barrier penetration using a 
Monte-Carlo code. The 3d equations of motion (PT|) were 
re-derived for an arbitary orientation of the lopsided dis- 
tortion. Starting from some randomly-chosen initial val- 
ues, an orbit was evolved in this fixed potential for a 
time tcoh- The orientation of the lopsided distortion was 
then randomized and the integration was continued in 
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X = log,o(1-e) 

FIG. 10: Angular momentum distributions from Monte-Carlo 
simulations in which the potential was re-oriented suddenly 
(top) or smoothly (bottom) eachfcoh- Red (rightmost): a — 2; 
black: 5 = 4; blue (leftmost): a = 8. Dashed lines show the 
predicted barrier location, Eq. (|50|l . 



the new potential, followed by another randomization of 
the potential etc. In addition to the parameters Ad, Am 
defined in Eqs. (gS]) and (gT]), this Monte-Carlo model 
has the additional parameter 

R = i^oUh = « 4a-3/2, (59) 



the dimensionless time between potential reorientations, 
where the final expression uses Eq. (|47bl) . The number 
of Schwarzschild precessional periods between potential 
re-orientations is ~ i?/27r(€)^. 

Monte-Carlo experiments were carried out for the fol- 
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FIG. 11: Time-averaged angular momentum distributions of 
stars in the Series III A'^-body integrations, in three intervals 
of semi-major axis. The distributions were computed using 
all stars with instantaneous 2 values in a range AlogjgS = 
±0.05 centered on the stated value, over the time interval 
< t < 2 X 10*^ yr. The solid (black) curves exclude stars that 
eventually become EMRIs; the dotted (blue) curves include 
these stars. Cross-hatched (grey) areas show the predicted 
location of the Schwarzschild barrier, Eq. H36|) . given the lower 
and upper limits on a. Hatched (blue) areas show the capture 
angular momentum for EMRIs, Eq. (|62|) . Solid rectangles 
show the angular momentum at the assumed capture radius 
around the MBH. 



15 



lowing sets of parameters: 

Am - 7, Ad = 3, R ^ lA (a w 2) 
Am = 30, Ad = 10, i? = 0.5 (5 w 4) 
Am = 120, Ad = 30, i? = 0.2 (5 « 8). 

For each choice of parameters, 1000 Monte-Carlo exper- 
iments with different initial seeds were carried out, and 
each experiment embodied 1000 re-orientations of the po- 
tential. 

Fig. llOl shows the resulting, time-averaged angular mo- 
mentum distributions. Also plotted there is the expected 
location of the Schwarzschild barrier, computed using 
Eq. ([5D|) . While the latter is by nature approximate. 
Fig. [To] reveals a tail toward low angular momenta rather 
than a sharp cut-off at any value of £; orbits sometimes 
reach values of £ that are ~ an order of magnitude smaller 
than the predicted ^sb- 

The angular momentum distributions in the TV-body 
integrations are shown in Fig. 1111 At all radii, there is a 
sharp cut-off in the distribution at some value of (1 — e). 
At large distances, a = 4 and 8, this cut-off lies close to 
£ = €gB, while at smaller radii the distribution extends 
beyond the expected barrier location (see also Fig.lH]). By 
comparison, while the angular momentum cut-off in the 
Monte-Carlo experiments is also quite sharp, it occurs at 
£ values that are somewhat lower than £sb for all values 
of a. 

The assumption that the potential changes suddenly 
every ~ tcoh is unrealistic. In reality, changes in the back- 
ground potential are due to the combined precession of 
individual orbits, which is a gradual process. One con- 
sequence is that adiabatic invariance will be respected 
for orbital actions whose conjugate angles are varying 
on time scales much shorter than tcoh- This is not the 
case if the potential changes instantaneously, as in the 
model just considered. For instance, if the period of 
Schwarzschild precession of an orbit is short compared 
with tcoh, its angular momentum will be nearly con- 
served. From Eq. ([5^ . this condition is satisfied for A^- 
body orbits near the barrier when 5 is sufficiently large; 
e.g. for a = 4, tcoh/tcB, ~ 60. For a < 1, this ratio 
is < 10, suggesting that adiabatic invariance will not be 
strictly enforced. This is a plausible explanation for the 
better success of the Monte-Carlo model at smaller radii. 

To test this idea, we carried out a second set of Monte- 
Carlo experiments in which changes in the background 
potential were continuous with respect to time. The total 
change in the orientation of the torquing potential after 
each tcoh was the same as in the first set of experiments, 
but now the vector describing the distortion was rotated 
at a fixed rate, along a great circle, from its initial to 
final orientations during each interval. 

Fig.lTUlshows the results. As expected, the angular mo- 
mentum distributions are now truncated more sharply at 
small values. In the case of a = 8, the distribution falls to 
zero at £ ~ £sb- As a is decreased, the distributions ex- 
tend progressively farther below the barrier, approaching 



more closely to the results of the first set of experiments. 
These distributions are quite consistent with those from 
the A^-body models. Fig. [TT] 

As noted above, the semi-empirical criterion p6p im- 
plies that there is a minimum value of a, 5, < Csb ~ 0.7 
using the adopted value of Csb, below which there is no 
barrier. A straightforward prediction is that stars with 
initial values of the semi- major axis below 0.7 mpc 
should be able to form EMRIs, at a rate that is unaf- 
fected by the arguments presented in this section. We 
show below that this is in fact the case. 

Nevertheless, a robust result of the work presented in 
this section is that resonant relaxation itself is ineffective 
at coaxing stars much past the Schwarzschild angular 
momentum barrier. Rather, these results imply that the 
barrier should be "hard," at least on time scales compa- 
rable with tpjpj or tcoh, or 10^—10^ yr in these simulations. 

C. Barrier penetration. II. 

As noted above (cf. Fig. [2]), one or two stars per Myr 
were captured by the MBH, on average, in the Series 
III integrations, most of them as EMRIs. Fig. [12] shows 
several examples. 

In Newtonian systems, classical, "two-body" (non- 
coherent) scattering is much less effective than resonant 
relaxation at changing stars' angular momenta when the 
motion is nearly Keplerian. But this is not necessarily the 
case for stars on orbits near or beyond the Schwarzschild 
barrier, where RR is effectively quenched by the rapid 
relativistic precession. In this section we consider the ex- 
tent to which classical, or non- resonant (NR), relaxation 
can explain the EMRI events in the A^-body integrations. 

The orbit-averaged NR relaxation time ^nr for stars of 
semi- major axis a in our model (a n cx 7'~^ density cusp 
around a MBH) is 

tNR ~ 4.6 Myr a^/^ 

with N the number of stars within 1 mpc (Appendix |B| . 

Suppose that NR were the only mechanism capable 
of changing stars' angular momentum leftward of the 
Schwarzschild barrier. The condition for capture onto 
a GW-dominated orbit would then be obtained by re- 
placing tnn by ^nr in Eq. 

tow = 2(1 - e)tNR. (61) 

We find from Eqs. (jlVB]), ^ and dM]) the critical value 
of a at which this condition is satisfied: 
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FIG. 12: Evolutionary tracks for a subset of the stars from 
Series III integrations that became EMRIs. Dotted (red) line 
is the Schwarzschild barrier, Eq. (|36|) . Dash-dotted (blue) line 
is Eq. (|62p . a prediction for the critical eccentricity at which 
GW energy loss dominates the evolution, assuming that the 
gravitational perturbations are dominated by non-resonant 
relaxation. Dashed (black) line shows the assumed capture 
radius. 

where N is the number of stars with a < 1 mpc. 

Eq. (15^ is plotted as the dot-dashed (blue) line in 
Fig. [T^ After crossing this line, stars can be seen to 
remain near to it for some time as their energy drops. 
This diagram suggests that Eq. ([62]) accurately speci- 
fies the region where GW energy loss and gravitational 
scattering are equally important, consistent with our as- 
sumption that NR relaxation is the dominant mechanism 
for angular momentum evolution in this part of the (a, e) 
diagram. 

Given a criterion for when a star enters the GW regime, 
we can then ask how often the barrier penetration de- 
scribed in the previous section would have resulted in 
EMRIs. 

Before doing so, we note two characteristic radii asso- 
ciated with flQw- When 

the Schwarzschild barrier lies to the left of the GW line. 
In the A^-body models considered here the critical value 
is a « 20, beyond Omax = 10. At and above such radii, 
the Schwarzschild barrier would not be an impediment 
to EMRI formation, and the EMRI rate (per interval 
of semi-major axis) would be similar to what was found 
above in the Series II integrations. This limit is proba- 
bly of only academic interest however, since in standard 



TABLE I: Mean times to EMRI formation in the Monte-Carlo 
experiments 



a 


ti (yr) 


t2 (yr) 


2 


3.1 X 10'' 




4 


2.8 X 10* 




8 


1.2 X 10* 


1.5 X 10* 



models of nuclei, almost all EMRIs would originate from 
orbits with a < 0.01 pc. 

At the other extreme in radius, the Schwarzschild bar- 
rier lies to the right of e = 0. In the A^-body models 
the intersection occurs at a = Csb ~ 0.7; in general, 
Eqs. ([35]) and (|62|) give for this condition 

/M 

a < 1.6 X lO^^pc f ^ j N-^^^. (64) 

Since the barrier does not exist at these radii, the dif- 
ferential capture rates would also be similar to what was 
observed in the Series III simulations. 

For values of the semi-major axis between these two 
extremes (0.7 < 5 < 10 in the A^-body models), the 
Schwarzschild barrier exists and lies to the right of the 
critical eccentricity for GW emission. EMRI formation 
at these radii requires a substantial degree of barrier pen- 
etration. 

We tabulated how often in the Monte-Carlo experi- 
ments from the previous section a star passed the GW 
boundary. Since not every experiment resulted in such an 
event, the mean event time, in each set of experiments, 
was computed using a formula from survival analysis [36j : 

7 I ^MC - A^e 

i=0 

where A^mc is the total number of experiments, A^e is 
the number of experiments in which the star satisfied 
the condition ([5^ at least once, ti is the time at which 
this first occurred, and T is the total elapsed time per 
experiment. 

The results are presented in Table 1, for the first (sharp 
changes; ti) and second (smooth changes sets of 
Monte-Carlo experiments. As expected, for large a, the 
two times are similar, since the barrier is no impediment. 
At smaller radii, the mean times are interestingly short 
only in the first set of experiments; in the second set, no 
events were observed. We note a certain "conspiracy" : at 
small a, the degree of barrier penetration is greater, but 
the GW line lies farther from the Schwazschild barrier. 

Next we consider the effectiveness of non-resonant re- 
laxation at penetrating the barrier. Several factors are 
relevant: 

1. Because RR is so rapid to the right of the barrier, 
the angular momentum distribution in this region should 
remain close to that associated with an isotropic phase- 
space density, i.e. N{t)M constant xiM. This con- 
trasts with the case [e.g. f2£|] where NR alone determines 
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the phase space density, leading to a logarithmic decrease 
in N with respect to £ near the loss-cone boundary. 

2. The angular momentum of a star near the barrier 
oscillates, at roughly the Schwarzschild frequency, with 
amplitude Si ^ £+ - i+ - Isb (Eq. Mil- To push 
a star past the barrier, a NR perturbation will require a 
finite amplitude A^nr ^ S£. 

3. Stars remain near the barrier only for a time ~ 
tcoh', after this, the direction of the background torque 
changes, and the star random-walks to larger angular 
momenta, as discussed in the previous section. 

The change in £ due to NR over an interval of time 
equal to tcoh is 



{A£) 



NR 



/ ^coh 
V^NR 



1/2 



(66) 



This change is large enough to move a star leftward of 
the barrier if 



NR. 



> 



^£4 



'£+-£^^2{e)^Ar>. (67) 
for which this con- 



Let fmax(a) be the largest value of j 
dition is satisfied. Writing 

£+ + l {e+ ^£-)^{£) + {ifAo (68) 
and eliminating {£) in Eqs. (|57)) and (|68p then gives 



1 



\l/2 
JNR 



(69) 



At every a, we expect stars with £^ < £max(a) to be 
scattered leftward of the barrier in a time tcoh- The frac- 
tion of stars at a with ^sB(a) < £ < ^max(a) is 



F{a) 



^max 



(a) 



(70) 



and the time scale for stars to be lost past the barrier is 
therefore 



tioss(a) 



1 dN 
N~dt 



i^(a)-4coh(a). 



(71) 



We evaulate (A£)j^j^ using each of the two choices for 
tcoh discussed above: the mass precession time, Eq. (|51ap . 
which gives 



(A£) 



NR.M 



4.4 X 10" 



(72) 



in the iV-body models; and the vector resonant relaxation 
time, Eq. ([5^ . for which 



{A£) 



NR.RRv 



5.3 X 10 



(73) 



Tables HI] and |TTI] give the computed values of F and 
-F'~^tcoh for the two choices of tcoh- Predicted loss rates 
are similar for high a values, and in both cases, NR is 
predicted to fail to breach the barrier when a < 2 — 3. 

Fig.HHplots histograms of the capture events in the Se- 
ries III integrations. As predicted, the number of events 



20 



15 



10 




-0.5 0.5 1 

log [a/mpc] 

FIG. 13: Distribution of semi-major axes for the capture 
events from Series II (top) and Series III (bottom). Red 
(unfilled) histogram shows the plunges; blue (cross-hatched) 
histogram shows the EMRIs; the total is indicated in black. 
In the upper panel the initial value of a is used; in the 
lower panel, the value of a during the final crossing of the 
Schwarzschild barrier was used. In both panels, the elapsed 
time is 2 X lO'' yr. To the left of the dashed vertical line in the 
lower panel, non-resonant relaxation is predicted to be inef- 
fective at pushing stars past the Schwarzschild barrier. To the 
left of the dash-dotted vertical line, the Schwarzschild barrier 
does not exist. 



falls sharply for a < 2 — 3 mpc. A few captures also occur 
from orbits with a < I, roughly the minimum value for 
which the barrier is present. 

We can also compare predicted and measured event 
rates. From Fig. [21 the mean capture rate at early times 
in the Series III integrations is ~ 1 — 2 x 10~^ yr~^- 
(Four out of 19 of the events were associated with or- 
bits below the Schwarzschild barrier, reducing the mean 
rate of barrier-crossing events slightly.) The number of 
stars initially with o>2 — 3is~ 35 — 40, and the loss 
times in Table 2 are roughly 1 x 10^ yr at these radii. 
The predicted event rate is therefore 3 — 4 x 10~^ yr~^ - 
in reasonable agreement with the measured values given 
the crudeness of the model. We note that our model can 
be expected to overestimate the capture rate since it ig- 
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TABLE II: NR loss rates: t^^h = 



5 t 



coh 



(yr 



2 1.2 

3 
4 



xl(F 
1.6 X 10* 



<-max <-SB 



(yr) 



4.9 X 10~^ 3.1 X 10* 

1.8 X 10* 2.1 X 10~^ 8.5 X lO** 

5 2.0 X 10* 2.5 X 10"^ 8.1 x 10® 

6 2.2 X 10* 2.5 X 10"^ 8.9 x 10® 
8 2.5 X 10* 2.2 X 10"^ 1.2 x lO'^ 
10 2.8 X 10* 1.9 X lO"'^ 1.5 X 10^ 



TABLE III: NR loss rates: tcoh = tRR v 



tcoh (yr 



F — ' 



(yr) 



2 2.6 X 10* 

3.9 X 10* 5.6 X 10^^ 7.2 x 10® 

5.2 X 10* 7.2 X 10"^ 7.4 x 10® 

6.5 X 10* 7.5 X 10"=* 9.0 x 10® 

. 7.8 X 10* 7.3 X 10"^ 1.1 X lO'^ 

8 1.4 X lO-""' 6.9 X 10"^ 1.5 x 10^ 

10 1.3 X 10^ 6.6 X 10"^ 2.0 x 10^ 



nores the possibility of a star returning to the right of 
the barrier after crossing it. 

We will present a more detailed calculation of the bar- 
rier penetration rate due to NR in a later paper (ssj . 



VI. DISCUSSION 

Here we discuss briefly how the key results from the A'^- 
body experiments can be extended to nuclear star clus- 
ters with more general properties. We treat this topic in 
more detail in Papers II and III [s^ [13]; in particular, 
we do not attempt here to derive absolute EMRI rate 
estimates for general clusters. 

We begin by collecting some of the important relations 
derived above and expressing them in more general form. 

Combining the parameter dependence of Eq. (|35|) with 
the empirical normalization of Eq. p6p . we find for the 
angular momentum that defines the Schwarzschild bar- 
rier: 



CsbV fr,V fM,V 1 



0.7 J \a J \ m J N 

here and below, N is the number of stars within radius a. N(a) cx a was not assumed in deriving this expression. 
However, that assumption was made in deriving Eq. ()62p . the condition that GW emission dominate stellar encounters. 
We can generalize that relation to a cluster with arbitrary density profile using the approximate scaling of the NR 
relaxation time: 



^NROc^j^ (75) 



MlP, 

N 

[e.g.[l^, together with the exact, orbit-averaged expression for tNR in the case n{r) cx r^^, Eq. ([BU)) . to write 



The condition for GW emission to dominate relaxation then becomes 

-1 / A f \ V5 / ^ ^-2/5 ^ ^ X -2/5 

(77) 



ej^^^^xiu y^^^j {iO^MqJ \WMqJ yW 

If physical capture by the MBH is assumed to occur when r < rcapt, ?'capt = ©''gj then the critical eccentricity for 
capture is 



We define Opiungc to be the value of a such that (l — e^)^^^^. = (l — e^)^^; at larger a, all stars plunge. We find 
that apiungc is defined implicitly by 

«« . 1,1 X io> (I)"'" f^y , i™. 
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FIG. 14: Illustrating the critical curves defined in the text for nuclear star clusters obeying density laws n(r) oc r"'', with 
various slopes and normalizations. Dashed (black) line: capture radius (vp = 8rg); dash-dotted (blue) line: radius at which GW 
emission dominates stellar perturbations, eq. (|77[) . The Schwarzschild barrier, eq. ((74]), is shown as the red line; it is solid where 
conditions allow EMRI formation. Below the horizontal line, non-resonant relaxation is expected to be inefficient at pushing 
stars past the barrier (eqs. ISOllSTllSSl ) M, — W^'Mq and m = WMq were assumed. 



The Schwarzschild barrier intersects the GW line when 



3/5 /r^ \ 2 / ,r \ 13/5 / \ -8/5 

120^ ^) ^) (80) 



and it intersects the capture line when 



(^)(jL)^,„o(^y(?v'(^)'( 



VmpcyVloV V0.7y \8 J \10^MqJ \10Me 



(81) 



One of these two relations defines the effective upper limit to the radial extent of the Schwarzschild barrier. Setting 
e = in Eq. (TM)) gives the lower radial limit: 



mpcj VlOV '^'^ V 0.7 j [iO^Mq) \10Mq) ' ^^^^ 

Another key parameter is the minimum value of a for which non-resonant relaxation is able to penetrate the 
Schwarzschild barrier (Sect. Vc). Combining Eqs. (|51ap . (pS)) . (I5D)) . and (|74p . we find that the critical value of 

a satisfies the implicit relation 

/ Q \ /C* \ / AI \ / ^ \ — 3/2 / jy \ — 1/2 

li)pe„e.ate"'H^j(l0^j (lO^j [w^) ' ^^'^ 

In deriving this expression, we have equated tcoh with ^m- 
We apply these expressions to nuclear star clusters obeying 

n(r) = nor-''' (84) 
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N = N{< a) = N^ih^"^ 
where d = a/mpc and A'^<i is the number of stars with a < 1 mpc. Combining Eqs. (j85p and 



, we find 



^penetrate ' 



140 



Csb\ ( M, 



0.7 J VlO^A^G 



5/2 



lOMo 



-3/2 



1 2/(5-7) 



N 



-1/2 
<1 



(85) 



(86) 



Fig. [14] plots the relations defined above for clusters 
of TO = 107^0 BHs around a M, = IO^Mq MBH. We 
chose likely values of 7 and A^<i following the discussion 
in Ref. [1^. As in the A^-body models, there is gener- 
ally a rather small range of a values from which EMRIs 
can form: small enough that GW emission can overcome 
stellar perturbations, but large enough that non-resonant 
relaxation can push stars past the Schwarzschild barrier. 
Interestingly, for sufficiently dense clusters, this range 
can go to zero, implying essentially no EMRIs; however 
it appears that the required densities are one or two or- 
ders of magnitude larger than expected for real galactic 
nuclei 

VII. CONCLUSIONS 

1. A^-body integrations have been used, for the first 
time, to directly simulate the long-term evolution of rel- 
ativistic clusters of compact stars around massive black 
holes (MBHs), both Schwarzschild and Kerr, and to com- 
pute the rate of extreme-mass-ratio inspirals (EMRIs). 

2. When relativistic terms are omitted from the equa- 
tions of motion, stars are scattered into the MBH at rates 
that are in good agreement with those expected from the 
theory of resonant relaxation (RR). 

3. Relativistic precession suppresses RR, leading to 
an effectively maximum value of the eccentricity at each 
value of the semi-major axis. This "Schwarzschild bar- 
rier" strongly inhibits EMRI formation, leading to cap- 
ture rates that are factors ^ 10 — 100 lower than in the 
non-relativistic case. 



4. We use an approximate Hamiltonian formulation of 
the perturbed equations of motion to explore two pos- 
sible mechanisms for barrier penetration: one related to 
resonant relaxation and the other to non-resonant relax- 
ation (NR). We show that NR is effective at penetrating 
the Schwarzschild barrier only for orbits with semi-major 
axes above a certain value, and this prediction is verified 
in the A'^-body integrations. Approximate expressions for 
the capture rate are derived and shown to be consistent 
with the rates observed in the simulations. 
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Appendix A: Hamiltonian model 

Here we use standard techniques [e.g. l26l| to derive the equations describing the rates of change of the Keplerian 
(osculating) elements of a star moving in the potential (1571) : 

CM / r \ 

^(r) = ^ + <i>p, $p^$,lnf — j -5'acos6l (Al) 

and including the time-averaged effects of Schwarzschild precession. 

We begin by transforming from Gartesian coordinates to Delaunay variables [e.g.Is^ which are action-angle variables 
in the Kepler problem. The Delaunay action variables are the radial action / = (GM.a)^/^, the angular momentum 
L, and the projection of L onto the z axis Lz- The conjugate angle variables are the mean anomaly w, the argument 
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of the periapse uj, and the longitude of the ascending node J7. In the Keplerian case, five of these are constants; the 
exception is w which increases linearly with time at a rate 

Vr {GM,f/I^. (A2) 

The Hamiltonian, averaged over w, is 



$p, (A3a) 



$p = ^^$p = i-^ dS(l-ecos£;)$p(r). (A3b) 

In the final term, E is the eccentric anomaly, where r = a(l — e cos E) and the eccentricity is e = -^/l — l? jP . After 
the averaging, l-i is independent of w, and / is conserved, as is the semi-major axis a. We are left with four variables 
and with $p as the effective Hamiltonian of the system. 

The orbit- averaged Hamiltonian describes slow, precessional dynamics. Superposed on the slow variations described 
by the averaged dynamics are fast oscillations, with frequencies ~ {GMt/a^)^/'^ and with fractional amplitudes 
(5 « a ($ — fGM,. If (5 <C 1, i.e. if M, ^ M,^, we can ignore these fast oscillations. 

After expressing the Cartesian coordinates in terms of the Delaunay variables, the results of the averaging are 

= $M + *Z5 + $GR, (A4a) 

$M = + F(e)] , (A4b) 

a 

C(a) = In (^^^ + 1 -ln2, F{e) = In (l + Vl - (?) -\J\~ e^, (A4c) 

= 5*06 sin i sin cli, (A4d) 

*GR = ■ (A4e) 

The averaged dipole potential, is expressed in terms of the orbital inclination i where cosi = Lz/L; z = for an 
orbit that is perpendicular to the major axis of the dipole. The last term, $gr, reproduces the orbit-averaged rate of 
Schwarzschild periapse advance, Eq. (jl6p . The longitude of the ascending node, VL, does not appear due to symmetry 
of the potential about the z-axis. 

In the limit e ^ 1, F{e) -> -f/2. 

We define a dimensionless time r = v^t where 

3GM. 

= Vr — 5 , (A5) 

the Schwarzschild precession frequency in the limit e — >■ 0. Dropping constant terms (including terms that depend 
only on semi-major axis a), the dimensionless Hamiltonian describing the perturbed motion becomes 

= = - il - e^y^'"^ + AMF{e) + ADesinisinuj, (A6) 

with ^M,^D defined in Eq. (I43p . The equations of motion, 

duj__dH_ m _ _dH_ dn _dH_ de^ _ _dH_ 

dr ^ d£ ' dr ^ ^ du ' Ot ^ d£,' dr ^ ^ on ^ ^ ' 

are given explicitly in Eqs. ()4ip . 



Appendix B: Non-resonant relaxation 



Here we summarize the orbit-averaged equations describing changes in angular momentum due to non-resonant 
relaxation (NR) and derive the angular- momentum diffusion coefficient for the A^-body models [e.g. Isol [39ll40j. 
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In terms of the binding energy per unit mass E = — f^/2 + ip{r) = GM,/2a, where ?/'(r) = GM,/r, and the 
normalized angular momentum R = IL\ = l"^ — 1 — e^, the Fokker-Planck equation describing diffusion in angular 
momentum due to NR is 



dN _1 d 
~dt ~ 2aR 



(Bl) 



where 

N{E, R)dEdR = N{a, e)dade (B2) 

is the number density of stars in (energy, angular momentum) space, and ((Ai?)^) is the diffusion coefficient in R, i.e. 
the sum, over a unit interval of time, of (Ai?)^ due to encounters. 

Taking the limit i? — >■ and averaging over one orbital period, this becomes 

^=^dRV'dR) 
where jl{E) is the orbit-averaged diffusion coefficient: 

(B4) 

J Vr 2R 

and the integral is over one full radial period. 'p{E) is precisely the orbit average of the inverse angular momentum 
relaxation time defined by Hopman & Alexander 41] and henceforth we write = ^nr- 
Let f{E) be the phase-space number density of stars; it is related to N{E) by 

f{E) = -±-(GM.)-'e'/'N{E) (B5) 

= foE^-^/^. (B6) 

The latter expression assumes n{r) cx r^^; the A'^-body models have 7 = 2. The local diffusion cofficient is expressible 
in terms of f{E) via 

Inn ^ (3/,/. - /3/. + 2/0) , (B7a) 

Io{E) = / fiE')dE\ (B7b) 



I„/,{E, r) = {2 [V^(r) - E]}-'^'^ f {2 [^(r) ~ E']}'''^ f{E')dE' (B7c) 

J E 

where InA w ln[A/,/(2TO)] is the Coulomb logarithm; in the A^-body models InA sa 9. 
The orbit averages are 

setting 7 = 2, the value in the A'-body models, we find 

ME) = ^foE-^GM,)\ (B9a) 
li/2{E) = ^(lnl6-2)/o£;-2(GM.f , (B9b) 
hME) = ^{n-l2ln2)foE-'{GM,f (B9c) 



and 



3/1/2-/3/2 + 2/0 = CJoE-^ {GM,y\ (BlOa) 

C2 = ^(12 In 2-1) (BlOb) 
w 1.18533 (BlOc) 
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so that 

HiE) = C2^^G'm'lnAf{E). (Bll) 

We note that Jl{E) cx f{E), a result that holds for arbitrary 7. 
Eqs. ([U IB2[ IBSp combine to give / in terms of a: 

/(«) = ^ (GM,)-'/' (B12) 
where Nq = r~'^N{< r) = a^'^N{< a). Then 

I 
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